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I.  INTRODUCTION 


Generally,  the  numerical  approach  to  finding  internal  wave 
eigenfunctions  involves  numerical  integration  of  the  differen¬ 
tial  equation  over  the  entire  ocean  depth.  For  a  deep  ocean, 
particularly  when  the  horizontal  wavelengths  of  interest  are 
much  shorter  than  the  ocean  depth,  numerical  integration  can  be 
confined  to  the  first  few  hundred  meters  of  the  ocean.  This 
simplification  arises  from  the  fact  that  for  a  wide  range  of 
possible  profiles  the  variation  in  depth  dependence  is  confined 
to  the  upper  few  hundred  meters  below  the  surface,  after  which 
the  profile  may  be  approximated  by  a  decaying  exponential.  More¬ 
over,  for  horizontal  wavelengths  much  shorter  than  the  ocean 
depth  the  exponentially  decaying  profile  model  may  be  extended 
to  infinity,  i.e.,  the  ocean  may  be  assumed  to  have  infinite 
depth.  Since  the  solutions  of  the  differential  equation  for 
a  continuously  decreasing  exponential  profile  are  known  in 
closed  form  (they  are  Bessel  functions  of  the  first  kind),  nu¬ 
merical  integration  of  the  differential  equation  needs  to  be 
carried  out  only  over  the  upper  few  hundred  meters  of  depth, 
within  which  an  essentially  arbitrary  Valsala  frequency  profile 
may  be  prescribed,  as  derived  from  experimental  data. 

The  numerical  procedure  to  be  presented  here  will  provide 
any  mode  function  and  dispersion  relation  over  any  wavelength 
range  for  an  arbitrary  Valsala  frequency  profile,  as  long  as 
it  decays  exponentially  at  great  depths.  The  practical  limita¬ 
tion,  of  course,  is  the  amount  of  computer  time  one  is  willing 
to  expend  for  results  covering  an  appropriate  dynamic  range. 

For  ordinary  ranges  of  physical  interest,  the  required  computing 
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time  appears  to  be  of  the  order  of  minutes  for  a  dispersion 
curve  and  seconds  for  a  mode  function. 

Dr.  P.  A.  Selwyn,  who  was  kind  enough  to  review  the  first 
draft  of  this  paper,  has  called  attention  to  the  fact  that 
there  already  exists  a  computer  program  (Refs.  1  and  2)  that 
performs  calculations  similar  to  those  undertaken  here,  but  for 
the  case  in  which  the  ocean  depth  is  finite.  One  would  expect 
that  it  might  have  been  possible  to  modify  the  earlier  program 
so  as  to  adapt  it  to  the  infinite  depth  case  as  well. 

Although,  in  a  sense,  the  question  of  that  possibility  was 
moot  by  the  time  this  information  had  been  disclosed,  neverthe¬ 
less,  the  option  still  exists  for  a  potential  user  who,  if  con¬ 
vinced  that  the  earlier  program  has  sufficient  advantages  over 
the  present  one,  might  be  willing  to  risk  whatever  numerical 
complications  may  be  involved  in  such  a  task.  Therefore,  it 
seems  proper  to  discuss  briefly  some  of  the  major  similarities 
and  differences  between  the  two  programs.  Unfortunately,  a 
direct  comparison  of  their  computational  accuracies  and  speeds 
is  not  feasible  at  present.* 

Both  programs  rely  on  the  same  general  method,  which  is  to 

start  with  the  appropriate  boundary  condition  at  the  deep  end 
••  ••  •• 

of  the  Vaisala  frequency  profile  and  numerically  integrate  the 
differential  equation  over  depth  for  selected  values  of  the 
wave  number  and  phase  velocity.  These  parameters  are  varied 
according  to  some  iteration  procedure  until  the  value  of  the 
solution  of  the  differential  equation  at  the  ocean  surface  is 
sufficiently  close  to  zero. 

The  earlier  program  solves  for  the  phase  velocity  in  terms 
of  the  wave  number.  Here  the  converse  is  true. 

*Tf  the  adaptation  of  the  earlier  program  were  based  on  the 
same  idea  used  here  to  account  for  the  boundary  condition  at 
infinity,  i.e.,  matching  with  the  Bessel  function  solution  for 
an  exponential  profile,  the  running  time  of  the  program  would 
be  increased  considerably. 
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The  earlier  program  uses  a  fourth-order  Runge-Kutta  method 
for  solving  the  differential  equation,  so  that  the  error  is  of 
fifth  order  in  the  incremental  step  size.  Here  the  Numerov- 
Manning-Mi liman  method  is  used,  and  the  error  is  of  sixth 
order  in  the  incremental  step  size.  This  gain  in  efficiency 
occurs  because  the  present  method  takes  advantage  of  the  fact 
that  the  differential  equation  does  not  contain  a  first  deriv¬ 
ative  term. 

The  earlier  program  uses  a  Newton-Rapheson  iteration 
scheme  to  find  the  eigenvalue,  whereas  here  the  regula  falsi 
scheme  is  used.  Generally,  the  Newton-Rapheson  iteration  con¬ 
verges  faster  but  may  fail  to  converge  if  the  trial  value  is 
not  sufficiently  close  to  the  true  value.  However,  once  star¬ 
ted  properly  the  regula  falsi  iteration  is  guaranteed  to  con¬ 
verge  .  * 

Other  differences  between  the  two  programs  primarily  have 
to  do  with  how  they  handle  the  problem  of  mode  jumping,  which 
occurs  in  the  calculation  when  dispersion  curves  for  succes¬ 
sive  modes  approach  each  other  too  closely.  The  details  are 
rather  involved  and  it  is  difficult  to  assess  the  relative 
merits  of  the  two  approaches  without  considerable  further 
study.  However,  it  may  be  worth  noting  that  the  nature  of  the 
regula  falsi  method  permits  the  use  of  certain  eigenvalue  prop¬ 
erties,  derived  from  the  classical  Sturm-Liouville  theory, 
that  aid  in  resolving  the  mode-jumping  difficulties  encountered 
by  the  program  presented  here. 

The  mathematical  statement  of  the  problem  to  be  solved  and 
some  properties  of  the  eigenfunctions  and  dispersion  relations 
needed  in  the  subsequent  development  are  summarized  in  Chapter 
II.  In  Chapter  III  the  general  framework  for  the  numerical 

•There  are  cases  in  which  the  convergence  is  extremely  slow, 
however.  When  this  happens  it  has  been  found  expedient  to 
switch  to  internal  halving,  which  is  less  efficient,  in  gen¬ 
eral,  but  which  converges  at  a  predictable  rate. 


solution  of  the  dispersion  curves  is  presented.  The  detailed 
description  of  the  algorithm  for  the  numerical  solution  of  the 
differential  equation  and  the  associated  eigenvalues  is  given 
in  Chapter  IV.  The  computer  implementation  of  this  algorithm 
INTMODE  is  described  in  Chapters  V  and  VI.  A  more  detailed 
description  is  given  in  the  Appendix. 

Some  examples  of  numerical  results  obtained  with  INTMODE 
are  presented  in  Chapter  VII. 


II.  DESCRIPTION  OF  THE  PROBLEM  AMD  SOME  GENERAL 
PROPERTIES  OF  INTEREST 
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Interest  here  is  confined  to  internal  waves  with  horizontal 
wavelengths  of  at  most  a  few  kilometers  so  that  the  effects  of 
the  inertial  frequency  do  not  enter.  Accordingly,  the  mode 
functions  are  determined  by  the  differential  equation 

ym(y)  =0,  (1) 

along  with  the  boundary  conditions 


d \p 


dy 


+  K 


N2  ( y ) 
2 
m 


-1 


(0)  =  0,  lim  (y)  =  0,  (2) 

m  y->  +  oo  m 

where  the  depth  y  is  measured  from  the  surface  y  =  0.  In  (1) 

K  is  a  given  wave  number,  is  the  mode  frequency,  and  N(y) 

ffi 

is  the  Vaisala  frequency  profile  associated  with  a  vertical 
thermocline  in  the  ocean. 

For  a  given  value  of  K  the  differential  equation  (1)  sub¬ 
ject  to  the  boundary  conditions  (2)  determines  uniquely  the 
infinite  set  of  mode  functions  ( y ) .  For  each  mode  a  disper¬ 

sion  relation  is  thereby  determined;  that  is,  each  is  a  well 
defined  function  of  K, 


over  the  interval 


n  =  ft  (K) ,  m  =  1,2,  .  .  .  , 
mm’  ’  ’ 


0  <  K  <  °° . 


(3) 
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j  3. X‘ G  'Jl n  G "V  —  _  ’ '  C  G  Ur  i"’ 
hat  vanish  at  X  =  0.  As  X  fceccr.es  -rfci- 
asymptoticaily  approaches  a  single  pcsi- 
to  the  maximum  value  of  N(.y),  the  Vaisala 


frequency.  As  K  approaches  zero  the  slope  of  each  of  the  dis¬ 
persion  relation  curves  defined  by  (3)  approaches  a  different 
finite  positive  value  l/um;  i.e.,  ^im  ^  =  ^im  ^  =  l/vm. 


The  dispersion  relations  (3)  can  be  expressed  parametric¬ 
ally  in  terms  of  a  quantity  y,  whicfci  approaches  um  as  K  approaches 
zero;  i.e.  , 


Km  "  Km^> 


ft  =  SL 


K  (y) 


m 


(4) 


where  the  set  of  functions  Km(y)  are  determined  by  the  differen¬ 
tial  equation 


d^ip 

— +  y2N2(y)iJ;  =  K2ip  ,  (5) 

dy2  K  w/ym  nrm* 

subject  to  the  boundary  conditions  (2).  Prom  this  point  of  view, 

for  each  value  of  y  that  is  chosen  above  the  minimum  value 

needed  for  the  existence  of  the  particular  modes  being  considered, 

the  boundary  value  problem  determines  a  set  of  positive  real 
2 

eigenvalues*  Km  (y),  the  square  root  of  which  provides  the  K  (y). 

According  to  the  general  theory  of  eigenfunctions  for  the 
linear  second  order  differential  equation  (Ref.  4),  if  N(y)  is 
piecewise  continuous,  bounded  and  approaches  zero  as  y  becomes 
infinite,  for  fixed  y  there  are  a  finite  set  of  real  positive 
eigenvalues  . *  As  y  increases,  the  number  of  these  eigenvalues 
also  increases,  but  for  y  smaller  than  some  critical  value  there 


_ _ _ _  _  o 

^Actually ,  the  convention  is  to  regard  the  quantities  -  K  as 
the  eigenvalues,  which  results  in  their  being  characterized  as 
negative  real. 
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r.ay  te  r.c  eigenvalues  at  all.'  It  is  alsc  kr.cwn  from  the  general 
theory  that  each  is  an  increasing  function  of  u  (cf.  Ref.  5, 
p.  357). 

It  is  assumed  that  the  Vaisala  frequency  profile  K(y) 
decays  exponentially  at  great  depths.  This  behavior  may  be 
characterized  more  conveniently  by  actually  requiring  that  after 
some  depth  yQ,  II ( y )  becomes  an  exponential  function;  i.e., 

_  Z 

N(y)  «  V  b,  y  *  yQ.  (6) 

As  observed  in  Ref.  1  the  solution  to  the  boundary  value  problem 
associated  with  (1)  and  (2)  for  an  exponential  profile  of  the 
form  (6)  is  proportional  to  a  function  <J> ( y )  given  by 

_Z 

4><y)  =  JK5  (MbN0e  b).  (7) 

Thus,  the  boundary  conditions  (2)  can  be  restated  as 

_Z 

MO)  =  0,  iKy)  =  A  JKb(ybN0e  b),  yayQ,  (8) 

where  A  is  a  constant  that  can  be  chosen  arbitrarily  or  to 
satisfy  some  normalization  condition. 

According  to  Ref.  3,  a  W.K.B.  solution  of  (1)  can  be  used 

1 1"1 

to  determine  the  approximate  m  J  mode  dispersion  relation;  I.e., 
within  the  W.K.B.  approximation  the  dispersion  relation  is 
determined  by  the  equation 

/  pdy  =  tt (m  -  -m) ,  (9) 

yN>a 


where 


P(y)  -  |  [N2 (y)  -  n2] 


1/2 


(10) 


III.  THE  NUMERICAL  PROCEDURE  TO  DETERMINE  DISPERSION  RELATIONS 


In  this  chapter  the  principles  underlying  the  numerical 
procedure  to  be  used  for  calculating  the  dispersion  relations 
will  be  outlined.  The  details  of  the  procedure,  itself,  will 
be  presented  further  or.. 


Since  K  is  an  increasing  function  of  y,  it  follows  from 
(9)  and  (10),  by  considering  the  limit  as  K  approaches  zero, 
that  the  minimum  permissible  value  y  for  the  m  mode  is  given 
approximately  by 


m 


(m  -  })tt 
~  ■  11 

J  N(y)dy 


(11) 


Since  (11)  is  only  an  approximate  formula,  in  order  to  guaran¬ 
tee  that  a  solution  Km  exists  for  the  eigenvalue  problem,  values 
of  y  chosen  to  calculate  the  dispersion  relation  for  a  given 
mode  should  be  somewhat  larger;  e.g.,  it  would  be  prudent  to 
confine  the  choice  of  y  to  values  such  that 


j,  +  (m-l)Tt 


jT 


N(y)dy 


(12) 


The  smallest  value  of  y  given  by  (12)  ought  to  be  large  enough 
to  guarantee  the  existence  of  a  real  Km  but,  ideally,  small 
enough  to  exclude  the  existence  of  any  higher  mode  eigenvalue 
Kn>  n  >  m.  It  will  be  found  that  this  is  generally  true  but 
that  there  are  some  noteworthy  exceptions. 
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For  any  choice  of  positive  K  and  y  it  is  possible  to  find 

a  solution  of  the  differential  equation  (5),  subject  to  the 

second  condition  in  (8)  v/hich  determines  the  necessary  initial 

values  to  be  imposed  for  that  purpose  at  the  point  yQ.  The 

differential  equation  can  be  integrated  numerically  from  yQ 

down  to  zero,  where  a  value  for  HO)  will  thus  be  acquired.  If 

the  maximum  magnitude  of  the  solution  over  the  interval  0  s  y  *  yQ 

is  then  a  function  4'(y,K),  determined  by  this  process, 

ms  x 

may  be  defined  by 


V(P,K) 


HO) 


max 


(13) 


As  defined,  ¥ ( p , K )  is  a  function  of  y  and  K  alone;  i.e.,  it  is 
independent  of  the  normalization  constant  A.* 

Because  there  are  cases  in  which  (12)  does  not  lead  to  a 
satisfactory  initial  value  for  y  (some  modes  are  skipped)  a 
slower  but  safer  procedure  than  relying  on  the  W.K.B.  approxi¬ 
mation  has  been  adopted  here.  The  differential  equation  (5) 
is  solved  numerically,  subject  to  the  second  condition  in  (8), 
with  K  set  equal  to  zero  and  y  set  equal  to  a  sequence  of  values 
— oo — - .  For  each  value  of  y  in  the  sequence  the  sign  of 

Jo  N(y)dy 

Ho)  Is  observed.  When  a  change  In  the  sign  occurs  the  cor¬ 
responding  value  of  y  is  used  as  a  trial  value,  and  it  is 
assumed  that  zero  bounds  K  from  below.  If  it  is  found  that 
modes  are  still  skipped,  the  y  increment  is  decreased  and  the 
procedure  repeated. 

2 

If  K  happens  to  be  an  eigenvalue,  then  *F(y,K)  vanishes. 
Thus,  the  problem  of  calculating  the  dispersion  relations  is 


•Because  it  automatically  relates  error  to  a  specified  dynamic 
range,  the  normalization  (13)  is  needed  for  stability  of  the 
numerical  process  used  to  solve  (14). 
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there  is  just  one  root,  and  for  y  in  the  interval 

U2  s  y  <  y^ 

there  are  exactly  two  roots,  etc. 

To  calculate  the  roots  of  (14),  interval  halving  or,  for 
more  rapid  convergence,  the  regula  falsi  (Ref.  6)  method  can 
be  used.  The  numerical  procedure  given  here  does,  in  fact, 
rely  upon  the  regula  falsi  method  to  solve  (14),  although  the 
technique  of  interval  halving  is  used  in  certain  circumstances, 
to  be  described,  in  order  to  reduce  computing  time. 

It  is  necessary  to  begin  with  two  trial  values  for  K(y), 
and  K2  >  K^,  such  that  and  4'(y,K2)  differ  in  sign, 

to  guarantee  that  the  root  K(y)  lies  between  and  K2;  i.e., 

K1  £  K(y)  s  K2  . 


Since  K(y)  must  be  positive,  initially,  the  trial  value  can 
be  zero.  For  the  initial  upper  bound  K2,  a  quantity  defined 
by 


K, 


N 


(15) 


where  N, 


max 


2  "max' 

is  the  largest  value  attained  by  N(y),  will  suffice. 


•The  fact  that  K~  as  defined  by  (15)  is  an  upper  bound  can  be 
seen  by  multiplying  (5)  by  and  integrating  from  0  to  <*>. 
Integration  of  the  derivative  term  by  parts  shows  that  that 
term  is  negative. 
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For  the  hither  modes  it  will  again  be  necessary  to  begin 
with  =  Q.  However,  may  be  set  equal  to  the  previously 
calculated  value  of  K  for  the  mode  one  step  down  at  the  same 
value  of  y.  That  is,  since  it  is  known  that 


Kn-l(u> 


Km<u), 


in  calculating  Km(y)  a  value  for  K2  given  by 


K2  <  Km_  q ( U )  (16) 

can  be  used. 

As  indicated,  the  trial  value  K?  should  be  slightly  less 

s  t 

than  Km_.(u)  to  avoid  accidentally  falling  back  onto  the  m-1 

mode  dispersion  curve  because  of  normal  errors  to  be  expected 

in  the  calculation.  A  test  should  be  included  here  to  guarantee 

that  the  choice  of  K~  is  not  too  much  less  than  K  ,(u):  the 

c  -  m- l  ~ 

function  y ( y , K )  must  change  sign  in  going  from  to  K2-* 

In  order  to  obtain  a  trial  value  satisfying  (16),  it  is 

necessary  to  have  an  estimate  of  Km_1(u)  that  is  known  to  be 

too  small.  Since  it  is  generally  not  the  case  that  K  .  would 

have  been  calculated  previously  for  exactly  the  value  of  y  now 

encountered  in  the  mode  m  calculations,  the  estimate  of  K  , (y) 

m— 1 

must  be  determined  by  interpolation,  e.g. ,  between  values 
Km-l(un)  and  Km-l(un  +  where 


However,  if  the  K  versus  y  curves  are  concave  upward  such  an 
interpolation  will  produce  an  estimate  that  is  too  large;  hence, 
the  desired  sign  change  in  '?(y,K)  would  not  occur.  On  the  other 

•While  the  regula  fatei  method  may  still  work  even  if  this 
requirement  is  not  met,  it  is  not  actually  guaranteed  to  con¬ 
verge  unless  the  sign  change  rule  is  imposed. 
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hand,  a  cruder,  one-sided,  interpolation  that  does  guarantee 

the  sign  change  can  be  used.  That  is,  instead  of  interpolating 

between  K  ,  (u  )  and  K  .  +  Ay)  the  trial  value  estimate 

m- 1  n  m- 1  n 

becomes 


K 


2 


n 


(17) 


From  the  fact  that  =  Km_1(y)/y  is  a  monotonic  increasing 

function*  it  can  be  readily  inferred  that  K2  defined  by  (17) 
will,  in  fact,  be  smaller  than  Km_^(y). 

Once  K  has  been  calculated  for  the  initial  choice  of  y  for 
a  given  mode,  the  value  of  y  is  increased  by  adding  a  small 
increment  Ay.  A  corresponding  increment  for  the  lower  bound 
trial  value  K,  can  be  obtained  from  an  interpolation  analogous 
to  (17). 

That  this  can  be  done  so  that  K1  continues  to  be  a  lower 
bound  can  be  seen  as  follows.  By  definition, 

K  =  yfi. 

Therefore , 


^  -  n  + 

dy 


dn 

dy" 


(18) 


Although  Q  increases  monotonically  with  y,  it  is  uniformly  bounded 
by  the  maximum  Vaisala  frequency;  hence,  the  second  term  on  the 
right  of  (18)  approaches  zero  as  y  becomes  arbitrarily  large. 

This  is  evident  in  view  of  the  fact  that,  since  the  derivative 
of  log  y  is  jj,  ~  approaches  zero  faster  than  i.  Then,  for 
large  enough  y,  according  to  (18), 


A  K  *+  A  p  • 


(19) 


•Ti1  it  were  not,  a  case  ot  anomolous  dispersion  would  be  implied 
since,  as  already  observed,  K  is  an  increasing  function  of  y. 
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Moreover,  because  Q  is  a  mcnotor.ic  function  of  u ,  the  estir.ate 
of  AK  given  by  (19)  is  always  too  small. 

Thus,  the  new  can  be  chosen  in  accordance  with  (19); 

i.e. , 

K 

K,  =  K  +  Ay  =  K  +  Ay  -±L  ,  (20) 

1  U  U  U  U 

where  is  the  previously  calculated  value  of  K  corresponding 
to  the  value  of  y  before  the  increment  Ay  is  added.  For  modes 
higher  than  the  first  (m  =  1)  the  use  of  (17)  to  obtain  the 
upper  bound  K 2  continues  each  time  y  is  incremented,  while  the 
lower  bound  is  obtained  from  (20).  For  the  first  mode,  how¬ 
ever,  (15)  is  the  only  estimate  immediately  available  for  the 
upper  bound  K2  as  y  is  incremented,  although  (20)  can  still  be 
used  to  estimate  the  lower  bound  K^. 


IV.  DETAILS  OF  THE  NUMERICAL  PROCEDURE 


A.  SOLUTION  OF  THE  DIFFERENTIAL  EQUATION 

The  Nunerov-Nanning-Millman  method  (Ref.  6,  pp.  204-205) 
is  particularly  convenient  for  solving  (5)  numerically  for  given 
values  of  y  and  K.  The  method  requires  two  starting  values;  for 
a  step  size  h,  Myo)  and  <4(yo+h)  must  be  furnished  initially. 
Then  the  differential  equation  can  be  integrated  by  means  of  a 
single  recursion  relation  that  involves  only  4>  and  its  second 
derivative,  which  is  obtained  from  \p  and  the  relationship  sup¬ 
plied  by  the  differential  equation,  itself. 

The  starting  values  of  ij;  are  obtained  by  recognizing  that 
at  yQ  and  yQ+h  the  profile  N(y)  is  an  exponential  function  of 

the  form  (6).  Thus,  in  accordance  with  (7),  at  these  points 

-L 

A  ^ 

ip(y)  can  be  set  equal  to  JKb(ybNQe  ). 

When  N(y)  is  prescribed  numerically  over  an  interval  (0,yo) 
the  resolution  of  N(y)  implies  a  limit  on  how  small  the  step 
size  h  may  be  taken.  Conversely,  a  natural  limitation  on  how 
large  h  may  be  is  the  requirement  that  it  be  small  compared  to 
the  minimum  wavelength  X  to  be  considered.  Since  the  wavelength 
is  given  by 


this  means  that  the  size  of  h  is  governed  by  the  largest  value 
to  be  considered  for  the  wave  number  K. 


15 


ate 


B.  SOLUTION  OF  THE  EIGENVALUE  EQUATION 

The  eigenvalues  K  that  determine  the  dispersion  relation 
for  each  mode  are  found  by  solving  (14)  over  an  appropriate 
range  of  values  for  y.  For  this  purpose  the  regula  falsi  method 
(Ref.  6,  pp.  4-5)  seems  most  effective. 

In  some  cases,  convergence  of  the  regula  falsi  method  is 
too  slow.  Therefore,  if  twenty  iterations  occur  without  satis¬ 
fying  the  prescribed  error  criterion  the  computer  program 
switches  to  interval  halving  with  an  error  criterion  applied 
to  K  rather  than  i(i. 

C.  SELECTING  INCREMENTS  OF  y 

In  accordance  with  (12),  the  increment  6y  used  to  obtain 
the  starting  value  of  y  in  going  from  the  dispersion  relation 
for  one  mode  to  that  for  the  next  is  normally  given  by  adding 
increments 


6y 


1 


N(y)dy 


(21) 


until  4i(y,o)  changes  sign.  At  the  start  of  the  mode,  is  then 
set  equal  to  zero. 


If  the  increment  Ay  along  a  single  mode  is  too  large,  a 
jump  to  the  next  mode  may  occur.  This  can  be  guarded  against 
by  calculating  ip(y  +  Ay,  K  )  which  in  that  case  would  have  a 
different  sign  than  i^(y,K^),  where  is  a  lower  bound  used  in 
calculating  Ky. 

Evidently,  as  a  practical  matter  Ay  must  not  be  too  large. 
It  is  also  true,  however,  that  Ay  must  not  be  too  small.  While, 
theoretically,  trial  values  are  chosen  so  as  to  guarantee  the 
necessary  sign  change  In  ^(yjK)  for  the  regula  falsi  method, 
in  practice  it  turns  out  that  when  Ay  is  sufficiently  small  the 
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sign  change  nay,  nevertheless,  fail  to  occur.  This  is  due  fc 
the  residual  calculation  error  in  the  K  that  corresponds  to  the 
value  of  p  before  it  is  incremented.  This  error  is  sufficient 
in  some  cases  to  overcome  the  theoretical  inequality  relied 
upon  in  the  derivation  of  the  rule  for  selecting  K2* 

A  compromise  rule  for  selecting  the  p  increment  is  to  let 
Ap  be  about  yq  ^p.  That  is,  a  reasonable  choice  that  seems 
adequate  in  practice  is  given  by 


Ap  = 


D.  ESTIMATING  THE  ERROR  IN  K 


2/  N(y )dy 

J0 


(22) 


The  test  used  to  determine  when  to  stop  the  regula  falsi 
iterations  in  calculating  K  is  the  condition 


|Y(p,K)  !<e.  (23) 

The  value  chosen  for  this  purpose  in  current  applications 
is  10-7,  which  is  intended  to  provide  at  least  a  60  dB  dynamic 
range  for  the  corresponding  mode  functions. 

Therefore,  the  error  in  K  is  not  given  directly;  however, 
it  can  be  estimated  by  linear  extrapolation.  If  Yn  is  the  value 
of  V ( p , K )  that  just  meets  the  test  (23)  and  Yn_^  is  the  value 
of  Y ( u , K)  in  the  iteration  just  before  that  one,  then  the  quan¬ 
tity 


K  -  K  , 
AK  _  n  n-1 


AY 


Y  —  Y  . 
n  n-1 


(2H) 


where  Kn  and  Kn_^  are  the  corresponding  estimates  of  K  in  the 
two  Iterations,  is  approximately  the  rate  of  change  of  K  with 


I 
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in  K  corre- 


r 


r 


respect  tc  a  change  in  1  ( y  ,  K ) .  Then  the  error  t ^ 
sponding  to  e  will  be  given  approximately  by 


e 


K 


(25) 


The  error  estimate  can  be  used  to  prevent  the  anomaly 
mentioned  earlier,  that  too  small  a  choice  of  Ay  can  result  in 
a  failure  to  obtain  a  sign  change  in  ^(UjK)  using  the  trial 
value  obtained  by  means  of  (20).  The  idea  is  to  make  sure 
that  the  error  in  the  calculated  value  of  K  is  always  negative, 
i.e.,  that  the  calculated  value  of  K  is  too  small.  This  can  be 
done  by  subtracting  after  the  iterations  for  K  are  completed. 


V.  COMPUTER  REALIZATION  OF  THE  ALGORITHM 


The  computer  program  DISPER  was  designed  to  calculate 
the  K  =  K  ( y )  relationship  using  the  numerical  techniques  des¬ 
cribed  earlier  in  this  paper.  This  program  will  write  the 
(K,y)  pairs  as  calculated  along  each  mode  to  disk  or  tape  and 
will  plot  a  graph  of  the  (K,Q)  curves,  referred  to  as  disper¬ 
sion  curves. 


A.  INPUTS 

The  inputs  to  the  program  are  of  two  types:  (1)  para¬ 
meters  read  in  under  a  NAMELIST  option,  and  (2)  data  points 
read  in  from  punched  data  cards. 

1 .  NAMELIST/PARAM/XO,  B,  STOPK,  STOPMU,  ND,  EPS 

XO  Real 

The  X-coordinate  of  the  last  data 
point  of  the  numerically  defined 
function  N(X).  For  the  STD  data 
XO  =  220.  m. 

B  Real 

Decay  constant 

For  the  STD  data  B  =  1300. 

STOPK  Real 

The  maximum  K  value  for  which  the 
user  wants  dispersion  curves. 

STOPMU  Integer 

The  number  of  dispersion  curves 
to  be  calculated. 

ND  Integer 

The  number  of  data  points  +1  to  be 
read  into  the  array  N. 
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EPS 


2 .  Data  Cards 
ITITLE 


N 


B.  OUTPUTS 

1 •  Printout 

a.  The  parameters  defined  by  the  NAMELIST  option  are 
listed  at  the  end  of  the  program  for  verification 
purposes . 

b.  The  value  of  /*N(y)dy  is  printed  next,  followed  by 
the  values  of  N. 

c.  At  the  end  of  the  calculations  for  each  mode  the 
number  of  (K,y)  pairs,  the  maximum  estimated  error 
for  K,  and  the  complete  list  of  (K,y)  pairs  for  that 
mode  are  printed. 

d.  Occasionally  the  error  in  K  cannot  be  estimated. 

When  that  occurs  a  message  indicating  this  fact  and 
the  current  values  of  K  and  y  are  printed. 

2 .  Disk  or  Tape 

The  values  of  the  NAMELIST/PARAM/B,  XO,  STOPK,  EPS,  ND, 

STOPMU  are  written  on  TAPE2  under  the  format  (^E22.7>  215//). 

Next  the  values  of  N  are  written  to  TAPE2  under  the  format 

(8E10.5) • 


Real 

The  error  criteria  imposed  on  the 
numerical  solutions  to  the  differ¬ 
ential  equations.  For  the  STD 
data  EPS  '=  l.E-7. 


Integer 

Ten  character  title  of  the  N(x) 
values . 

Real  array 

Dimensioned  500,  read  in  on  punched 
cards  under  the  format  (8F10.5). 

N  contains  the  equispaced  data 
points  that  numerically  define  the 
function  N(X).  The  points  are 
spaced  a  distance  of  XO/(ND-2) 
apart . 
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At  the  end  of  each  rr.cde  the  number  cf  (K,y)  pairs  calculated 
is  written  tc  TAPE2,  format  ( // I 5 ) .  The  (K,y)  pairs  are  then 
written  to  TAPE2,  format  (2E22.7). 

TAPE2  may  be  defined  as  a  permanent  file  by  using  a  cata¬ 
log  control  card,  or  it  may  be  defined  as  a  magnetic  tape  by 
using  a  label  control  card. 

3.  PLOT 

A  10"  by  10"  graph  consisting  of  the  STOPMU  different 
curves  of  the  (K,fi)  pairs  is  plotted  at  the  end  cf  the  program. 

C.  EXTERNAL  REFERENCES 

DISPER  references  several  external  subroutines  that  must 
be  provided  by  the  user  through  control  cards  that  attach  the 
appropriate  permanent  files. 

The  necessary  routines  are  listed  below  under  the  name  of 
the  permanent  file  on  which  they  reside. 

1 .  INTMODE 

INITIAL  Reads  the  data  values  of  N(X), 

calculates  the  integral 

N (x )  dx 

Calculates  "best"  estimate  of  K 
given  y. 

Numerically  solves  the  differential 
equation  using  the  Numerov-Manning- 
Millman  method.  Called  by  GUESS. 

Function  to  find  MU0  for  each  mode. 

Writes  the  (K,y)  pairs  to  TAPE2 . 

Sets  up  the  calls  to  the  CalComp 
plotting  routines. 


GUESS 

DIFF 

XMU0 

OUTPUTK 

PLOTER 
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PSCALE 


Scales  the  axes  to  the  calculated 
data.  Called  by  PLOTER. 


ERRPRO 


2.  I  DALI B 


PLOTS 

PLOT 

LINE 

DAXIS 

SYMBOL 

NUMBER 

ABRTJOB 


3.  BESSEL 


JBESS 

JAIRY 

GAMLN 


Processes  detected  errors  through 
a  call  to  ABRTJOB.  Is  called  by 
all  of  the  routines  on  this 
permanent  file. 


CalComp  plotting  routines  called 
by  PLOTER. 


Error  processor  that  generates 
TRACEBACK,  prints  error  messages, 
and  terminates  the  job.  Called 
by  ERRPRO. 


Routines  to  calculate  the  BESSEL 
functions.  Acquired  from  the 
Argonne  National  Laboratory. 


D.  ERROR  MESSAGES 

We  have  attempted  to  anticipate  some  of  the  errors  a  user 
might  encounter  when  using  DISPER  under  very  general  conditions. 
If  one  of  these  errors  is  detected  by  the  program  the  error 
processor  ERRPRO  is  called.  ERRPRO  does  three  things,  (1) 
prints  a  brief  message  describing  the  error,  (2)  indicates  in 
which  routine  the  error  occurred,  and  (3)  terminates  the  job 
without  a  dump. 

A  table  of  error  messages  generated  by  DISPER  and  possible 
corrective  actions  that  might  resolve  the  problem  is  presented 
below. 
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TABLE  1.  ERROR  MESSAGES  AND  POSSIBLE 


CORRECTIVE 

ACTIONS 

Message 

Si gni f i cance 

Action 

Issued  By 

END-OF-FILE, UNITS 

NO-1  is  greater  than  the  number 
of  data  points  provided  for  N. 

Decrease  ND  or  provide  more  data. 

initial 

INDEFINITE  OPERAND 
9/9 

Solution  to  differential  equa¬ 
tion  is  inconsistent. 

Loo*  for  coding  errors  in  the 
routine  DIFF. 

DIFF 

DIVISION  BY  ZERO 

Algorithm  for  solving  the  differ¬ 
ential  equation  has  broken  down. 

Data  points  may  be  too  far  apart. 
Introduce  more  data,  pernaps 
through  interpolation. 

DIFF 

FIRST  MODE  CANNOT 

BE  CONSTRUCTED 

Cannot  find  MU0  for  this  mode. 

Ree xami ne  data. 

XMU0 

TOO  MANY  POINTS 

More  than  500  (K,u)  pairs  are 
need  to  construct  this  mode. 

Decrease  STOP*,  or  revise  program 
by  re dimensioning  SAVEK,  SAVE MU, 

T E M P K  ,  and  TEMPMU. 

DISPER 

K1  and  K0  ON  SAME 
SIDE  OF  CURVE 

Estimates  for  K  are  not  upper 
and  lower  bounds. 

GUESS 

(1)  EPS  too  severe. 

(2)  Extrapolating  initial 
estimate  of  Kl  beyond 

K  values  of  the  preced¬ 
ing  mode. 

(1)  Reduce  EPS. 

(2)  Reduce  STOPMU  or 
increase  STOPK. 

E.  DETAILS  OF  JOB  EXECUTION 

The  following  is  a  sample  card  deck  for  executing  DISPER. 


1. 


|PU»  T n t V T | . 

LaHEL 'PLnT»PE,"*pT^®-t*NO«aEPLOT *Osin»X»«Wt vSu*0t 
ATlACH(REe«>FL.f«H*i  .In*pnj 
AT  f  ACM  ( fOAl  TP.ll»3r<5> 

AT  T  AC^n  I  *Op  t  A' •  10*Pn*MR*l  * 

MAP(OFP, 

L0«0(0l5PBt»l> 

LnScT ll rp»TPAL iM^PCSeEL/I^TfOoF' 
tE*tCUTE. 

78 

89 

fp ApAH  . nDs)  13  »5T0PK«,  I  *Fj»i«1  ,F-.7f 

piilsf 


,"i  o5ooono[ni('So*»(M*<roio5f’"0*q,<MO«o('oo'' 
,  "I  050000P  fll09l»P«M)O4Ol<H''"0'0,0l0?0$00" 

#fllo5ooflno.«^nsonflno*o10S''n(l•'0.91oB,'^,^0', 
4"l o^onnnn  flio^ooflrO“oio?nnO*o,VioBo«oo'' 
,f»lo50ono<l*n1  O^onnflO^OlO^nnO^o.yi  0«fl«no« 


.(M0«0nu**0."l»*90«l006 

S|85JJ88l8;.1i8t888R8;8i858SSR8;2jf}§38fiS8 

nin5noo5O.A'09P(»000.(M050BU«0,”ln9oOnnj 
fll  (l5(*flO;0.fl1(l5UBOOO.OlO*0««nO."lB5noOBQ 
nl  i>5n00,0.Oi05o«0n<1.0i0?0nOM>,',,ln900nn0 
r»i  n9«oo;o,nio5unOofl.Oio?flno«o,<,lflSoooon 
noA0T6?i^ 
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This  will  produce  a  listing  cf  the  ( K , y )  pairs  for  each 
mode  and  a  plot  of  the  dispersion  curves.  It  will  not  create 
a  tape  or  permanent  file  of  the  (K,y)  pairs. 

If  the  user  wishes  a  permanent  record  of  the  (K,y)  pairs 
additional  control  cards  must  be  included. 

For  a  permanent  file  the  deck  would  contain  two  additional 
cards.  The  sample  deck  below  is  an  example. 


IOU.To.mTi.  nPAPes.3,375* 
flF'JUEST  <TaDf?,*MF5 

LABEL  (Pi OTAP€,«,RTF6'.C,N0«oePl.0T,0BL0,xB«v,wSM.0, 

ATTACH(PE«*FL.MM*i . InaPD) 

ATTACMJIOAI  TP.IOSCG) 

ATTACH (nr sop IM. fQsPp.MRal  1 
MAp(OEF) 

LnAO(OISPPTN) 

LnSETIltpatnALlH/PFSoEL/rf'iTMooP1 

EXtCUTE. 

Catalog ( T«oF?.uiSPcUovE.inaPn) 

\ 
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If  the  (K,y)  pairs  are  to  be  written  to  a  magnetic  tape 
the  job  stream  might  look  like  the  example  shown  below: 
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Note  the  change  in  the  job  card  as  well  as  the  additional 
control  cards. 

Should  both  a  magnetic  tape  and  a  permanent  file  be  desired 
the  job  stream  would  look  like  the  next  example. 
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VI.  PROGRAM  MODE 


The  computer  program  MODE  was  designed  to  calculate 
and  plot  the  normalized  mode  functions  4'(y,K).  Given  (1)  the 
dispersion  curves  created  by  DISPER,  (2)  a  set  of  consecutive 
mode  numbers,  and  (3)  a  value  for  K,  this  program  will  use  the 
dispersion  curves  and  linear  interpolation  to  find  the  corre- 
ponding  y  values.  It  will  then  calculate  and  plot  the  mode 
function  for  each  mode  number. 


A.  INPUTS 

The  inputs  to  this  program  are  of  two  types:  (1)  the 
outputs  of  DISPER,  and  (2)  a  NAMELIST  option. 

1.  TAPE5 

TAPE  5  is  defined  to  be  the  disk  file  or  magnetic  tape 
produced  by  DISPER. 

B , X0 , STOPK , EPS , ND ,STOPMU  The  first  record  on  TAPE5 

is  the  defining  parameters 
used  by  DISPER  to  create 
the  dispersion  curves. 

They  are  read  in  under  the 
format  (4E22.7,  215//). 

For  definitions  see  inputs 
to  DISPER. 

N  Empirical  data,  format 

(8F10.5).  See  Inputs  to 
DISPER.  There  will  be 
ND-1  values  of  N. 

NPTS  Integer 

The  number  of  (K,y)  pairs 
for  the  current  mode. 
Format  (//I5). 
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XI,  MU1 


Real 

The  K  and  MU  values  of 
each  mode,  (2E22.7). 


STOPX  Real 

The  depth  to  which  the 
mode  function  is  to  be 
calculated.  Must  be 
greater  than  or  equal  to 
X0. 

IFIRST  Integer 

The  first  mode  to  be 
calculated. 

LAST  Integer 

The  last  mode  to  be  cal¬ 
culated.  All  modes  bet¬ 
ween  IFIRST  and  LAST  are 
calculated. 


B.  OUTPUTS 
1 .  Printouts 


a.  The  parameters  defined  by  TAPE5  and  the  NAMELIST 
option  are  listed  at  the  end  of  the  program  for 
verification  purposes. 

b.  N  Is  listed. 

c.  Mode  number  Is  printed  followed  by  a  list  of  PSI 
values  for  that  mode.  Format  is  (4E22.7). 

2 .  Plots 

10"  by  10"  graphs  of  the  (X,PSI)  values  will  be  plotted, 
one  plot  for  each  mode. 


C.  EXTERNAL  REFERENCES 

MODE  references  several  external  subroutines  that  must 
be  provided  by  the  user  through  control  cards  that  attach  the 
appropriate  permanent  files. 

The  necessary  routines  are  listed  below  under  the  name 
of  the  permanent  file  on  which  they  reside. 
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1.  I  NTMODE 


NITIAL 


DIFF 


Reads  the  data  values  of  N'X ' 
calculates  the  ir.terral 


:: ( x )  ax 


Numerically  solves  the  differen 
tial  equation  using  the  Numerov 
Manning-Mi liman  method. 


PSCALE 


ERRPRO 


INTNP 


FUNCT2 


FUNCT4 


PLOTMOD 


2.  I  DALI B 


PLOTS 

PLOT 

LINE 

DAXIS 

SYMBOL 

NUMBER 


Scales  the  axes  to  the  calcu¬ 
lated  data.  Called  by  PLOTMODE 

Processes  detected  errors 
through  a  call  to  ABRTJOB.  Is 
called  by  all  of  the  routines 
on  this  permanent  file. 

Calculates  the  integral 


I  N2  (x  )  tf2  (x  )  dx 

V 

Calculates  X  *  J(v,X)**2. 

Called  by  INTNP. 

Calculates  the  alternative 
asymptotic  approximation  for 
X  *  J ( v ,X ) **2 .  Called  by  INTNP. 

Sets  up  the  calls  to  the  Cal- 
Comp  plotting  routines. 


CalComp  plotting  routines 
called  by  PLOTMOD. 
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A3RTJ0B 


GAUSS 


3.  BESSEL 


JBESS 

JAIRY 

GAMLN 


D.  ERROR  MESSAGES 

A  table  of  error  messages  generated  by  MODE  and  possible 
corrective  actions  that  might  resolve  the  problem  is  presented 
below. 

TABLE  2.  ERROR  MESSAGES  GENERATED  BY  MODE 
AND  POSSIBLE  CORRECTIVE  ACTIONS 


Error  processor  that  generates 
TRACEEACK,  prints  error  messages, 
and  terminates  the  job.  Called 
by  ERRPRO. 

Numerical  integrating  routine 
using  Gaussian  quadrature.  Called 
by  INTNP. 


Routines  to  calculate  the  BESSEL 
functions.  Acquired  from  the 
Argonne  National  Laboratory. 


Messige 

SI gn 1 f 1 cance 

Action 

Issued  By 

More  than  500  pal rs 
were  needed  for 
this  node. 

DISPER  was  altered  to  permit  more 
than  500  pairs  to  be  calculated. 

Mike  similar  changes  in  MODE. 

MODE 

The  maximum  K  value 
on  TAPE5  Is  less 
than  the  K  of  In¬ 
terest. 

The  dispersion  curves  were  not 
calculated  to  this  value  of 

Rerun  DISPER  with  STOP*  greater 
than  K,  or  reduce  value  of  K. 

MODE 

There  are  not 
enough  modes  on 

TAPE  5 . 

STOPMU  Is  less  thin  LAST. 

Reduce  LAST  to  less  than 

STOPHU,  or  rerun  DISPER  with 

STOPHU  greater  than  LAST. 

MODE 

Indefinite  operand 

9/f 

Solution  to  differential  equa¬ 
tion  Is  Inconsistent. 

Look  for  coding  errors  In  the 
routine  DIFF. 

DIFF 

Division  by  Zero 

Algorithm  for  solving  the  differ¬ 
ential  equation  has  broken  down. 

Data  points  may  be  too  far  apart. 
Introduce  more  data,  perhaps 
through  interpolation. 

DIFF 

E.  DETAILS  OF  JOB  EXECUTION 

The  following  is  a  sample  card  deck  for  executing  MODE, 
when  the  dispersion  curves  are  on  a  permanent  file. 
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IPl). Ml  |  •  unl^F-i  t 

aTTaC"'  I  MUuF  * 

FTNt  l=i'OOE»L=u> 

MAP  <01- f) 

ATTaC- I  lNlMOut  •  iO  =  *-n  > 

ATT„c"ioeisei..Ti,*p,,i 

ATTACr<  i UAL  lb  «  I usC»> 

LAbtL  <  PLO  f  A?t  •  A  *H1  ’iij*L=Nob_,t.f'LUl  *0sLO*AsSV*VSN=u) 

LuStT  <l.lesINTrtnUE/„ES‘>Ei./  u,ALi0> 

LOU. 


* INPUT  IF IHST*1 ,L«sT*2 ,i\  =  .02,Sl UPA=^SU.k 


Next  is  a  sample  card  deck  for  executing  MODE  when  the 
dispersion  curves  are  on  a  magnetic  tape.  Mote  the  VSN  number 
should  be  the  one  assigned  at  the  time  DISPER  executed. 


I  HU  «  M  I  ^ 

aTTaL"  IMUUF  •  l.jsPU) 

Fil'd  1  =MOot  »L*"  • 

HAP  lUt  K  » 

*TTaCm  <  INl  -’OuF  «  1  U*K<> I 
ATTftcnfotP'iEL.  fu*pul 
ATTACH  I IU«l lb* IU*C«> 

LAbtL  <  TAPtSfH.NOHImi>»LsUl.->t'ClJH',t  *  VbN«bbbbl 
LAbtL  <PLU  TaPf  «  «*K  I  m,.«  |  »Nu5rtEP(_uT  fO*LO*A*5V  »  VSN*° ) 
LoStT ( L lo* INTHOof/nt SStL/ IDALIb) 

00. 


*  Input  1HBsT  =  1  «  LasT  *2  *n*.i>2*  SI  OP  AsdSO.  5 
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VII.  EXAMPLES  OF  DISPERSION  RELATIONS  AND  EIGENFUNCTIONS 

OBTAINED  WITH  INTMODE 

The  two  Vaisala  frequency  profiles  considered  in  the  sample 
calculations  employing  INTMCDE  are  shown  in  Fig.  1.  The  profile 
labeled  "exponentially  stratified  ocean"  corresponds  to  a  deep 
ocean  without  a  thermocline  and  K(y)  =  .00528  exp  -  y/1300 
radian/sec,  where  y  is  in  meters.  This  profile  is  identical  to 
the  one  used  by  Garrett  and  Munk  (Ref.  7).  For  the  exponential 
Vaisala  frequency  profile,  the  mode  functions  are  Bessel  func¬ 
tions.  Consequently,  results  of  INTMODE  for  this  profile  can 
be  compared  with  results  based  on  analytical  formulae,  thus 
providing  a  check  on  the  accuracy  of  the  numerical  technique. 

The  sharp  thermocline,  labeled  "STD  data  set",  is  taken  from 
(Ref.  8)  and  is  based  on  measured  towed  thermistor  chain  data 
in  the  tropical  Pacific  Ocean.  The  data  extends  to  a  depth  of 
220  meters;  at  greater  depths  an  exponential  profile  with  a 
decay  constant  of  1300  meters  is  assumed. 


The  dispersion  curves  for  the  first  25  modes,  correspond¬ 
ing  to  the  STD  data  set,  are  plotted  in  Fig.  2.  The  angular 
frequency  is  in  radians/sec.  For  the  STD  data  set,  plots  of 
the  first  four  internal  wave  modes  are  shown  in  Fig.  3  for 
K  =  .01  radians/meter  (X  ~  628m)  and  in  Fig.  4  for  K  =  .02 
radians/meter  (X  =;  328m).  The  mode  functions  are  all  normalized 
in  accordance  with 


/ 


^2(y)  N2(y)dy  =  l 


Since  N(y)  decays  exponentially  with  depth,  this  normalization 
constraint  leads  to  a  progressive  increase  of  the  mode  maximum 
with  mode  number  and  depth,  a  feature  corroborated  by  the  plots 
in  Figs.  3  and  4. 
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INTMODE  is  capable  of  yielding  mode  functions  of  any  order 
(subject  to  course  to  the  resolution  of  the  data  with  respect 
to  depth.  Figure  5  shows  a  plot  of  the  25th  mode. 

The  dispersion  relations  for  the  exponential  profile  are 
shown  in  Fig.  6;  the  first  four  mode  functions  are  shown  in 
Fig.  7  and  Fig.  8,  for  K  =  .01  and  K  =  .02  radians/meter, 
respectively. 

Examples  of  dispersion  relations  for  other  profiles  are 
shown  in  Figs.  9  and  10.  The  corresponding  profiles  are,  re¬ 
spectively,  those  referred  to  in  Fig.  12  as  STD  data  and  NRL 
data.  An  additional  example  of  dispersion  relations,  corre¬ 
sponding  to  a  pulse  shaped  profile,  is  shown  in  Fig.  11.  This 
profile  is  of  the  form 


N(y)  =  0; 

0  <  y  <  65. 45m, 

N(y)  =  .0105  rps 

65.45m  <  y  <  111m 

N(y)  =  0; 

y  >  111m. 

0.0250 
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MODE  NUMBER  t 


3. 


The  first  four  internal  wave  modes  for  the  STD 
data  set  (Fig.  1)  (X  cs  628  meters). 
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FIGURE  4.  The  first  four  internal  wave  modes  for  the  STD 
data  set  (Fig.  1)  (X  =*  314  meters). 
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FIGURE  9.  Internal  wave  dispersion  relations 
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FIGURE  10.  Internal  wave  dispersion  relations 
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FIGURE  11.  Internal  wave  dispersion  relations 
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GI„EN  A  F,,mC|TON  n,.. 

NUMEOICAI  1.V  from  *  .  0  to 
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(*<.*2  -  MO**2  v  N,X,*«p,  •  Hsi 
“ITh  PaRaPfTfRS  MU  ggi0  K  That  Ta^e 

on  various  rfal  values 


IT  IS  KNOWN  that  in  The  PEOIOm  Whehf 
x  .GF«  *n  THAT  ONE  SOLUTION  OF  THE 
DIFFERENTIAL  EOUaTIon  is  The  bE*SfL 


EFN?Ti0w  *JTw  ivf1k*  *„*  p  ANO  arum 
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xn  pEaL 

tHf  x.COOpniN*TE  OF  T^k 
l_ AcT  POINT  OF  THF  NllMEHT- 
C  Al  DaTa 

B  BEaL 

THf  DfCAV  constant 

NO  TNtEGPB 

THc  NUMBED  of  numeRtcaLl* 
DEfINfd  POINTS  ♦  1 

StOPK  BEaL 

thf  maximum  r  value  to  Rt 
CAICUiaTF.it  FOB  EaCh  HOOF 


STOPmU 


TNTEfiFP 

THf  number  of  modes  To  Rt 
CALCui  AtFO 

PE  aL 

fBoOB  CMITeMIA  IMPOBEf  OH 

thf  Solution  to  the  diffe¬ 
rential  FOUATJ0N 


The  last  IMPiiTS  to  the  pborram  aHf 
JHE  NyMRFpTCAL  vAL()Fp  THAT  DEFlKt 
N  ( X )  FBoh  n  TO  Vfl.  TMET  ape  Read 

fpom  caods.,  t«e  v»lues  must  bf  su 
Punched  aS  Tn  rf  head  unde0  'he 


Punched  aI’Td  rfEheSoU^n 
FobmaT i »fi n .s) 

N  bEaL  aPhav 

DImENfIUN  «oo 


: 

i 

L 


COMMON/ INPUTS/R.no,XS,Ep' ,RTPPmU, SToPk,  STOP* 

COMMON  PST I  S(j  Q) •N(tOj) *»'0*DFLlAX*NT, ITITLF 
Common  /Out/  Sa''FK(5oo>,  SaVfmU(50D) 
niMpNslON  temPk  I  son) ,TF“PMu (koO) 
pfal  N*  NS.  MU.  Kn.  HMAV »  Mllp.  K 1 
TmTeGEr  StoPmU 

n  aMel iST/P aR AM/o .NO. yq  *  *Topmij  »  STOP* » EpS 
WOOF  ■  1 
FoRMAX  ■  0. 

PT  m  ACOS l-j.» 

TfBo  ■  ft 


pcu^Pam  nlSPFR 


m 


1  ?o 


175 


135 


1*0 


145 


150 


155 


100 


105 


170 


C 

C 

C 

C 

c 

c 


•»ooo 


c 

c 

c 

c 

c 

c 

c 


10 

c 

c 


c 

c 

c 

c 


c 

c 

c 


?oo 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


PpAO  PARAMETERS 

0 FAc  param 
point  PaRam 

Pp An  nUMEOICAL  »AL"ES  F«R  N  ANO  INTpGOATt  N 

CALL  INIT|AL(SIIM1 ) 

POINT  lOOn*  SUMi 
PPRmaTI*  SUm1*,«*,f2p.7> 

POINT  lOO,  (N(It .T«1 «MO» 


WOITE  HAMfLIST  »n(J  N  TO  TAPP? 


CALL  OUTPnTK  M, 

FTNfl  MAXIMUM  V At  HE  FAR  m 

KMAXaN(l) 
oo  To  **2t»Nn 
KMAX»AmAXT(K“AX.NIt) ) 

CALCULATE  “UP  Ecu  THf  FtRST  branch 

nv  a  . l/SuMl 

muO  ■  *muo (O.tnMt 

CALCULATE  TNlTlAl  fSTIMaTE  FOR  THF  JNpReMfnT  IN  mU 
ALCnG  a  Ml  I p  branch 

ThElT  ■  ,K  /  SU“1 
xoElT  ■  TnELT 

START  The  PROCFAB  oE  APPROXIMATING  k 

Kfl  s  P? 

Mtl»MUO 
KCLNT  ■  0 

continue 

CALCULATE  the  UOPEP  FSTtMATp  OF  K,*l 
t<l  a  MlJ  *  KMAX 


CALL  SURRnUTINF  GUFSS  TO  l“OPOVF  The  FSTlMaTt  K 
IS  uppFp  fSTtmate.  k:  n  lo*<fR  fsTtm*If.  the 
fMPROVtO  f5TTHa7f  tS  RetUHMfc  JN  K», 


CALL  GUESS  |KC,K»  .Mil) 


CHECK  EOR  EPROP  cOOEb  RcTUhnfO  IN  NT 
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PROGRAM  r>  |  SPFo 


175 

lflO 

1«5 

190 

1,5 

ZOO 

Zo5 

210 

315 

ZZO 

ZZ5 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


£ 

c 

£ 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

£ 

£ 


c 

c 

c 


M  ■  0  T^nirATES  The  aLGORJTM  eOl'Ln  WOT  PpO 

D'lff  A  K  POM  MtiO  ON  THE  PINS  I  «RaNcH 
Kl  1,  mOt  A»OvE  THE  CuRvE. 


IP  (NT  ;EQ.  9 )  C  At  t  FRnPRo  (5> 


the  na*IM'im  eSTtmaTEo  E“RUh  poh  k  rs 
RETURNED  qV  0UF«*  IN  THP  X40TABLF  *o 


FPSk  a  KQ 

FflfiHA X  ■  AMAX1 (ppSk*fRP«A*» 


increment  kount.  ip  pqo  (ho,  k )  pairs  ma»f 


PFEn  CALCULATED 
E&HOR  MESSAGE  f 


«4ptBlJilErit  A  FT 


[OPK  NOT 


mm  §5 


GfAF»AlEO.  ANO  THF  Job  1WOBTEU.  EITHER  Bfe- 
nUCF  ST*PK  OR  TmpBfA«E  THfc  OlHpNSION  OF 
5AVFK (boot ,5»yF“tt(50S.)TFHPK(5oO) .TpMPHU (Sun) , 
Then  chaNoe  the  test  on  kount 


KOLnT  a  KOUNT  ,  1 

rp(K0*JNT  ,pe.  5««*  C  Ai  L  pRMPoo«IFR> 

S4Vfk (KOUNT)  »  pi 
*A VF^U  IKOnNf )  a  My 

CALCULATE  the  L«“ER  K  EsT1m»TE  FOR  THp  Kt«l 
HU  VAL^E 
increment  MU 

Kn  .  Kl  *  TOELT  •  Kl  /  mo 
If ( Kq  tot.  sTOPa>(jo  to  ino 

V(J  a  MU  ♦  TpptT 

96  FrBMAT<lX,T3;Z(«*tFZS.5>> 

GO  TO  ZOO 


HP  ARE  FINISHED  "ITH  THIS  HPANCH 

300  CONTINUE 
NPTrbkOUNT 

cURRpMT  <M'|.K>  PAT oS  To  USE  I*  ESTIHaUNG 
K 1  ALONG  the  NE*T  MOoE 

no  900  IaT, KOUNT 
TfNPK(1)»5AVFK(t< 

900  TpVpHU*I)a5AVEH"»I1 
CALL  OutpmTK  <KO,,NT) 


PRINT  IHE  maximum  fRpOR  FSTtmaTEO  FOR  K  Non 
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1 1  fir  i«  ifiaifciUi 


PRUt.RaW  nlSPFn 


230 

235 

240 

245 

250 

255 

260 

2*5 

270 

275 

250 

2*5 


c  rue  muoc 

c 

point  95,  mooE,  fpfk 

On  foRuaT<«  poTiMarro  Maximum  fPHoR  In  5  FO*  mooE  •»  14.  •  r<  •* 
1F?2.7> 

C 

C  BFGIN  1H£  NEXT  qOANCH 

MOCF  »  MOnf  *  1 
IF  (MODE  .rtf.  STopmii)  nO  TO  4*0 
1000  CONTINUE 
►T  a  0 
C 

woo. a  XMUfi  (M(iO,ny> 

KOan,. 

>  =  UU° 

FOUNT  *  0 
TrELT  ■  XnELT 
LrtCUNT  ■  „ 
pOj  CONTINUE 
C 

c  FIRST  APPROXIMATION  OF  a  FOB  THIS  BRANCH 

C  RFC-IN5  ST  FINOTWfl  THf  K  FOH  THE  LAST  BRANCH 

c  .  . 

00  ?0  I“1,NPTS 

TFITEMpHU(!)  .(3f.  mU)GO  TO  ,0 
20  C0KTIN(jE 

T  «  NPf S  *  i 

30  5|  CPE  *  TFMPKtI-j)  /  TE“PH"n-l> 

C 

c  calculate  the  upper  fsttmatf  of  k,*i 

c 

c 

*1  a  SLOPf  #  Mu 
C 

C  CaLc  SUBROUTINE  MJeSs  to  improve  The  FSTimnTfc  K 

C  5 ]  IS  UPPpR  ESTIMATE.  KS  15  LOato  fsTTMaIf*  Tmf 

C  IwPpQVEo  FSTiMATF  tS  PetUknfp  in  K1, 

C 

CALL  GUESS (K0fK, ^MU, 

c 

C  CHEc5  FOR  ERROR  rOOES  RfTunved  IN  NT 

C 

C 

c  1-T  .  <9  A  hRaNoH  me  ANTICIPATED  OTU  NOT 

C  OorUR  wERf.  A  result  OP  OltH  tsll- 

C  motion  PRoctnnPt  for  mUo 

c 
c 

T» (NT  .EG.  9) GC  TO  1"00 

C 

221  CONTINUE 
C 
C 

c 

c  The  MAXIMUM  FSTiuaTEO  FPPOH  foh  k  TS 

C  PfTiiRnEO  RT  rtUF««  tN  THf  vapIaHLE  *6 

C 
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4 


PBo(jB4m  MSPFo 


2<9o 

2<j5 

300 

3fl5 

310 

3*5 

3?o 

3?? 

330 

315 

3*0 


C 

BOSK  a  Ko 

FflfiMpX  «  «m4*i (rPSK>rRPu«X) 

C 

r 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 


c 

c 

c 

c 

c 


INCREMENT  kOHNT.  ip  sf)o  ,*ti,  p,  Ptl«$  H*vp 

«fEn  calc"Iateo  *Kr»  bTOok  rot  beached  an 

POBOB  HE5*AfiE  T<  PoImTeo#  t  pT*  fPOofl  tp 
flFNEB»f£D.  AMI  Twp  Jpifl  urUmteu*  EITHER  Pfc- 
Pf.'Cf  PJOP<  00  I»<rRFAB£  rbt  nt«€N5ION  Of 

«4VfK  <  boo »  .«*VPM1I(S0S»1  TF*PK(S00>  *TpHopulp«/»>  , 

kpLnT^OUmT  •  1 

IF<K0uNT  .06.  5«ft)P»lL  PPHppO |  If  R) 

SAVpK  (*OU»lT)  ■  V) 
bavfpu(koiint>  «  MU 

CALCULATE  THE  Lo*ER  *  F«TIm»TE  POP  TUP  Ntxl 
Ml.'  VALUE 
INCREMENT  HU 

K«  j  Kl  *  TPELT  •  K 1  /  “II 
TP<K<j  *6T.  STOP* i  ftO  TO  3oi» 

“I!  ■  HU  «  TOPLT 
00  TO  201 


*oo  continue 

L«u  a  SaVfmu ( komhT t 


pNOflLE  2 
bF*1N0  2 


C*LL  PLOTeR<nPTp,,,km*xi 


position  the  pup  to  tmp  bfoinninr  of 
The  PINST  more 

pPA0(?,95,«,*0,«TCPK tEP«,Nnt«TOPHU 
95  P0BM*T(*ET5.T*?ra//) 

PRINT  Pa«*H 
TpToP  «  Nn 

PpAn<?*9*>  (Nd<.  T  a  1.  I  ATOP) 

9*  FoRm jT  <  0F1 0 .5 ) 

or  ftoOO  J  a  It  arooMip 


c  NO*  TO  CALCl'LATp  0“EoA  FON  p*CH  HOOP  aNr 

C  P(.CT  0*£Ga  VS  k 
C 

C  I'p£  ApHAT  F4VE"m  TO  mOLo  Thf  uheGa  VA|  UfS 

c 

pFA0,2,  90) NPTa 
00  5o6»  I  a  i t  «PTS 
pfao(?*  1"0)  S««iFK(I>»  mii 
eAVfMUU)  a  pAyP*  ( T )  /  »U 

9?  foBm»T  I  //', l«) 

100  FoBm4T<?E22,71 
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PfluuR/lP  MSPFR 


101  FORMAT  (IX, 

Fnoii  fp^TINUe 

34S  r»u.  pi-0TfR(nPT«.2.  F«a») 

.006  rrMiNUF 

(*4Lt  PI-Otfo(NPt«»3«EpBm»*) 
r»LL  ENOP|  OT 
C 

3<50  C  jrB  is  CO«PLfTF 

C 

STOP 

F*'0 


smbOcut  !»,£ 


o.uESS 


I 


5 


10 


IS 


20 


?5 


1» 


15 


*n 


*5 


50 


55 


•n£CK  Ol'ESS 

Cy5S8N^iN5ll?^§?Nn?Js«^n^STOPMtUSToPK*  STUPA 

COMMON  PSj  (KflO)  .N(S0«)  *mO«ofi  1  **  t  NT  »  I T I  TLF 
0F4U  N0,N,M„,K1 .»8,5» 


c 

c 

c 

f 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 


8 

c 

c 


TF«R  «  2 

ThE  haMnhm  eST  tmjtepi  fobub  fun  k  t5 

BFTijBNEO  BY  fiUFSS  TN  THf  VaPTAMLF  Kg 

TF(Kg  ?Efl.  6.  )  FPSK  ,  S. 

SOLVE  OIFf  EQUAttON  uITu  hii.  Kg 

IF  K Q  CLOSE  ENOhbk.  OETmHN  with  KJ  a  Kg 


stTfmpI. To  f 

close*  To  *6 


!N? 


«N  UPPE*  BOUND  THAT  Is 
uTS  IS  more  FFFICIEHT  FOP 


Those  CaSfS  WHF»F  ko  15  A  CLOSER  ESU- 

n»Tf  THAN  K|. 


koLnT  ■  0 


K 1  ■  Kfl 


EPfJBO  to  1« 


Nr  EsT1m*tE  OP  tH«  EDp0o  IN  K  P0S5IRLE 


POINT  «8,  Kg ,  mm 

9ft  FCPmaTIo  no  ERRHO  fSTIM.TE  TN  K  POSSIbLF  MtPl**  /  •  K0 
1  5Xt  •  MU  ■  f2p,7» 

Kfl  a  ERSK 

PFTuBN 
in  continue 


TfkomX  a  aflS(PSTfjJ  f  PslMAX) 

Tf^PK  ■  Kf 
TFSf  «  PST (1) 

STEP  1  *  <K1  -  K«t  /  i^0« 

Kl  a  K» 

nr  POP  I  a  1 .  1" 

Kl  a  Kl  •  STFPl 

COLL  nlFF(M(J,KT.PSlN»Xt 

Tpfiu  a  ARS ( PS  I  1 T  T  /  DSImA*» 

TF(TEph  .be,  ep,,  no  To  ?0 

TfBmi  ■  AaS ( (Kl  -  TEmPK>  /  fTEWM  -  TEmPM**>  •  FPs 
IF  (TFHm1  .(IT.  FP5K )  FPS*  »  TfHMT 
Ko  ,  FRS* 

KT  a  Kl  -  TCOHl 
pfTmRN 
20  continue 

TfKpk  ■  K? 

TpKPMX  ■  TFPM 


A-ll 


■  *.  F2?.T* 


t 


phpoouiInE  *uE*5 


6ft 


65 


7ft 


75 


0® 


«5 


50 


95 


100 


105 


110 


20? 


C 

C 

C 


99 


200 


210 


22ft 

226 


236 

?36 

2*6 


c 

c 

c 


285 


Tp  <p$T  < 1 )  •  TE*t  .lT.  O.tun  TO  ?P1 
oftN  jiNuE 


pbRoR,  NO  5 1  ON  omAnGp 

POINT  99tMi),K0,a,.TEsT.i»«I»f» 

FORm|T11X,8f?0,»5| 

TP(KO  .NE,  6.)0«  TO  ?50 

ft- 7  ■  9 
OETmAN 


np*  USE.  The  OEpiM  *  F»L5T  hfthou 


continue 


*PL  NT  "0 
NC^ECO  ■  6 

nphfci  a  p 

K?*Ko  *  Te5T  •  f  K I «KS >  /  (TpeT  » 


KPLnT  ■  KOUNT  ♦  1 
IE(K0unT  .ST.  2"l00  tO  2*0 


OdVE  ■  PST (1) 

CALL  OlFF  (MIJ.K?.PS!HaX) 
tfpm  ■  RBgtPgipt  /  Psi“»x» 

TP  (TEhN  .fip.  E°5>  Go  To  2i S 

TpfiMl  ■  APS  < (K?  -  TEmPK)  /»TE«m  > 

TF(tE"m1  .Gt.  F»«Ki  rP**  a  tFh**1 

KO  a  EPSK 

Ki  a  K2  -  TEPNl 

RETURN 

CONTINUE 

TF*PK  ■  K? 


Tf*PPX  a  IfRa 

TF  (TEST  •‘psIIiM  2i0*?2«»«»?0 


CONTINUE 
CONTINUE 
KO  a  K2 

TpST  a  PSi ( 1 ) 
Pel  ( I  I  ■  e A VP 
GO  jc  200 
CONTINUE 
CONTINUE 
K»  a  K2 

o.o  TO  200 
CONjiNuE 
PAVE  ■  PST  (11 


INTERVAL  HALVING 


KOLnT  a  0 
CONTINUE 

koOnT  *  KOHNT  •  i 
TpIkOUNT  .rT,  ^*180  TO  -igO 


si»m 


TeMpMXI  ) 
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5IIO»rullN£  AUESS 


120 


1?5 


no 


135 


1*0 


1*5 


ISO 


155 


160 


165 


C 

C 

c 


Tr^PF  *  K? 

TrfPM*  ■  TfO« 

K)  a  (K]  .  K0)/», 

CALL  OlFF<MU,  K».  RS!MA»> 

TfRm  a  A0S (PS  I ( ) * /05  TMA*  > 

IF<Tf5m  •Lt«  Ep«'30  t0  ,0° 
tF«mU  A0S(K?-Tfm?i<>  /K? 

IF  (  TEpM  1  «LT.  l.r-6)Go  To  3ofi 
!F<Te5T  *  PST  1 1  *  .IT.  O.tGO  TO  PRO 

Kn  AND  K2  ABE  0*i  ThE  SA«F  RTnES 

K n  ,  k2 
TEST  ■  PSt  < 1 » 
r,r  TO  285 
28?  CONTINUE 

Fa  aNO  K2  ABF  Oh  OpPoSMT  5T0ES 

K1  m  K2 
SAWf  a  PS f  1 1 > 

00  TO  285 
296  CONTINUE 


c 

c 

c 

c 

c 

c 

C 

c 


FRRpR  PROCESSING 


AT  .  99 


IkioJcATES  OVER  50  ITERATIONS 
HAVE  OrCUoEU.  RUT  NO  SATISFACTORY 

FsTJmATE  oF  *  RtACHEP 


AT  a  99 
RETURN 

256  rrATINUE 


AO  SIGW  CHANGE  has  OrCURED.  EITHER 
ko  IS  TOO  IAPGF  oR  Ki  I*  loo  SHALL. 

THET  lIE  ON  TMF  same  SInF  Or  The 
oijRVE 

C»LL  ErRPbO(8) 

30P  rrKtlNUE 

aaRmaL  RETURN 

TFRM,  .  ARS,,K?  -  TEmPK^TfRH  .  TFMPmx,,  a  fcPS 

TF<TeB*1  . 0T •  fosaIEoSK  ■  TrRHj 

«n  a  EpSK 

FI  ,  k2  .  tFrHI 

RrlllRN 

F*  0 
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» 


I 


"fcCK  Wfl CE 


A  one  6*00  PhorRam 

• 

PROGRAMMER 

• 

P*  i.  nPaP^» 

• 

Im« T l TUTE  pom  npFFNSE  ANALYSIS* 

• 

A?l  TNflTnN*  Va> 

• 

HaooH  2s,  '<)«« 

• 

mCo£  C4LC1|E*tfS  a>iO  «P*Pue  mode 

• 

FUNCTIONS 

fop  otSpeR'IO*  CURVES 

• 

cbfateo  BY 

OamPaNtOM  PmorPam 

• 

OIFPES. 

• 

PROBLEM 

• 

6  T  WFN 

• 

(Ti 

THf  DTFFFPfNTIAL  pOHATinw 

• 

D»«?  tPsi>  /  ntx)  ••  9  ■ 

• 

(Ki*?  .  Mu**?  •  •  PSI 

• 

• 

NO  •  fXP  <- <X-Xfl| ^81 

• 

» 

Pf>H  x  .GE.  xO 

• 

(?) 

N(X)« 

# 

mIImfPTcaL  da*a 

• 

FOR  0  .LE.  X  .LF*  XO 

• 

(?) 

BOHNnApv  COmOItION 

• 

OSflO'  ■  * 

• 

(4) 

K 

(«) 

M(J 

• 

FTwo  THf  S"LUTtPN  TO 

• 

Tmf  OIFfeRfMIaI  EQUATION 

• 

c 

C  TKPllTS 


c 

NO 

InTfpEP 

NIJMbFB  OF  UATA  POInTR 

c 

R 

REA| 

OFCaV  con«)AnT 

c 

X„ 

REai. 

Last  x  ruoMDTN#TE  of 

c 

MIMfRIlaLi  T  nErlMEn 

c 

n#T»  k(x» 

c 

FPP 

REai 

FRPoP  LBITfcRiA  FOP 

r 

solution  of  oIffe- 

c 

OFNTIaL  tOUATIoN 

c 

K 

PEAI 

The  V  a  Li  it  of  K  UF 

c 

TNTfPFST 

c 

RTOPX 

PEAI 

MOW  Fah  01 1 1  THf  x  - 

c 

ax I c  MOne  functions 

c 

ape  To  hc  lalciiLaTfp 

c 

N 

PEAI 

apRaY  Imat  HOLnS 

c 

the  numfktcal  oxta 

c 

nt«ENSTONrn  Apn 

c 

irlRST 

tNTr^EM 

FIRfT  MOUF  OF  tNTE- 

c 

WFST 

c 

Last  ImTE 

1  aSt  MUnt  OF  ImTE- 

c 

"fSt 

A-l^ 


AO 


A  A 


70 


75 


50 


«5 


5(1 


55 


1  no 


ln5 


110 


c 

c 

f 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c 


£ 

c 

c 

£ 


c 

c 


c 

c 

c 

c 

c 

c 


C'llPUT 


o*»  Pfcm  aupa*  Th«T  hnLnS  THf 

“ODr  FUNCTION 


o»Mf KiStONrn  Ann 


f(*»«MCN/ INPUT#  /“.MO,  *A,F''5,«TrPMilt$TnP«f  STUPA 

rrAMCN  PST<5ofl>  f'o'OFl  TaX.NT.ITITLF 

raFAL  K*  Mil*  K‘l  *0,  Ui  «jM0  .  wUl.  MO 

r»TFGfH  ST0BwU 

NA*FLlST'TNPn t /a,  NOt  *:«  FP«»  IFIBgT.  I  AFl.  K«  plOPx 
PT  a  ALOS/-1.) 

TFfia5 


pfao  THE  toEnTtfvinG  Parana  Tpw5 
F0Cm£  IaPfs 


OFAn  (5 ,9A)  5«*«,5T0PK,rPS,nr .5T0PMII 


OfAQ  T^PUT 
K.STOP** IrtB5T.(  A?T 

PFAn  i^put 


PPInT  NaMflI5T/TmPI|T/  FnP  VrPlFI- 
CaTtON*  T mFN  PPtwT  N|X) 

point  Input 

IFjK.gJ.  5T0PK,  tai.L  EBoPHn  ,  tem*i  v 
rrtLAS1  .cT.  StadpII)  CAI  L  FPPHHOl  IFP**) 

call  initial  to  bead  tmt  vai  i>ts 

OF  N<X>  Off  T*Or  5 

oft  N(X)  paTa*  awo  CaLC"L*Tf  delta* 

call  InITial (SUM1 , 
ppint  loi . (n( n . ?ai tun) 

POINT  9g 

7«ToP  a  No  »  1 

PFLTAMU  a  PT  /  «(1M1 
if<t'fimst  .f<i«  nao  ro 


ATVAKCt  PdTNTEP  TO  Tme 

FTFoT  m0(1f  OF  TNTFPE5T.  !HO<1 


TOoP  a  IFIP5T  .  } 
no  «00  I  a  ] •  T*TQP 
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PPu<tP»H  •‘Opp 


ns 


l?o 


l?S 


no 


ns 


1*0 


*45 


isn 


1SS 


ISO 


US 


170 


neap <s*99i mpts 
or  soo  J  c  i ,  mots 
rpap  (S*  10o»  nUMWv.niiMuY 
SOO  rroTlNUg 
SOS  rrKTINuE 

rail  °<-OTs<U0*  «n) 

no  CIO  1  ■  IFIRst.  LaST 
C 
C 

c  sob  each  pf  thf  «pnEc 

c 


OFAP <S.99)NPTS 
K  ft  a  0  • 

MtlO  a  (1  «  .p )  •  ofLTaMii 

rn  52n  J  a  1,  MO 

c 

C  °FAn  E*CK  (K,  PMFQA)  PATO  To  FIND  K  INTf*VmI 

c 

nFAn(S,10O)Kl,Mnj 
TF<K1  *6T.  K ) 60  TO  5*0 

KO  a  <1 
«U0  a  Mill 
52o  CONTINUE 
510  CONTINUE 

c 

C  PC  LINEAR  TNTERoolaTtON 
C 

Sl.CPE  a  (Si  -  f  f MU >  -  MIIQ, 

VTNT  *  Kj  «  Stoop  •  uU i 
M'J  a  (K  -  YINT)  /  SLOPE 
PRINT  99 
pqINT  101  ,HlJtK 
C 

C  Solve  OIFrERENTtAL  EpUatIOn 

C 

POLL  nlFF-Md,  k.  PSI«AX, 
raLL  IntNP(k»  mii,  sUm-j) 

AwAXn  c  S^THa*  '  SORT  (S'iL*3t 
r,  ”  a  St0bx  /  HF|  T&x  *  ’ 

C  FRBpR  CONoTTIQm 

C  chec*  that  nijm  wot  Greater  than  thf  dt*ensipn  pf 


c 

c 

c 


S6r 


c 

c 

c 


TF(NUM  .GT,  son, CALL  ER°PHO,ItW) 

SILL  I*  The  rest  of  ns  I 

MT  ,  nUm  _  Np 
T  T  NPX  *  K  *  R 
AOG  a  HU  •  p  *  wo 

CALI  ja^SS(l?Nn»,  1«  aRp  •  CXH ( « (L  ♦  1)  •  UFLTaX 
PSI (Nn  ♦  L) » 

CONtINuE 

moBmalUE 


B)  . 
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PCu(iP*W  •<()DF 


ITS 


lfln 


1«S 


iso 


1S5 


nr  sTo  L  ■  1 1  K'hm 

ps  I  (L '  *  PS  I  (L )  /  SOPTteilMT, 

STn  roNfiNuE 
C 
C 

C  nr  PLOTTING  ANn  nUTPtiTT TNG  MFHE 

c 

POINT  98,  I 

POINT  1 0 1 • (PSI ( »T ) •  T I  -  1«  MJM) 
r*LL  Pl-OTMOn  (M.'ui  j  ,K.K  /  mu,  a*AXN  ) 
T«T»BT  «  J  ♦  I 
Cr 

C  itvaNcE  PniNTEP  10  NfXT  muof 

C 

nr  580  J  _  tsTadt|  NpTS 
PpAn<5»  i*o>nUMuv,nUMMY 
SfiO  rrNTlNU£ 

MO  rrNfiNUE 

CALL  E^OPLOT 
STCp 

101  PrSMAT ( 1 X ,  4F2? . 7 ) 

ion  FrRM4T<2E??,7> 

99  FrPMAT(//f TS) 

9q  FrBMAT<lHj,«  MOiaf  nUuBE15 
9 ft  FnPM4T <4e1 S.T,PTS//) 

FKT 
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sngioul  I K  F  r>  1  FF 


1 


5 


in 


1*5 


30 


25 


10 


15 


40 


45 


SO 


•n£CK  ''IFF 

ei'HoOulINF  nf  FF  ,Mt, ,  K  ,  Pe  T  M»  „  ) 

ff'VMON/  INOMTF/fl.tJO,  X;,f-»S*^T0PmU«5ToPk*  STUB* 
ri'KuON  PST<5fl0).*(50n)  <  mO.(5c|  TaX*NT»ITITLF 

'JFAL  NU  ,N  ,M(I  ,K 
C 

c 

c  This  is  Tup  MJVcoflv  _  MaNntmb  -  million  fFtwurj 

C  rrP  SOLVInr,  SECOND  On[)|ff5  LTk'FAM  DIFFERENT  T*l 

c  roCATlONS 

c 

F(Xj,K*»2  -  PU •«?  •  *##5 
POST  (  X  *  Y)  a  F  (  X  1  •  Y 
wrtss«NO 

TSTOP  ■  No 

c 

C  OFFIKE  FIRST  Two  VaLlipS  OF  B«1 

f 

TrPP\s^#fl 

C4UL  JtiESs  (TEPOI #1 *temo2*f*p I-0ELTaX  /  bl*PSI(l|) 
CALL  J“ESc  <tEmo»  ,1  ,tEMo?,PcT  t2>> 
TfFm»0ELTaX**2/i9, 
p«ImAX  «  aRs(Pst<2)) 

c 

C  SOLVE  The  EOUATtoN  8»CKuAHn« 

C 

no  ioo  lBi,l?Tno 
»kLm«?.*Pfi ( t”1 '-P5I ( I-9> 

XK'LM*XNUM»nPei  <►,  (Mn-f  *P5?  (1-2)1  •TER*' 
„M.MaXNUM4nP«!l  (M,Nn_T*1,  ,PlT  (  1_1  )  )  .TERP.lP. 
orA0N»l .-TpRM*F »n(nd-I*9>  * 

PS  I  (  I )  “XNllM/ptNOM 
TFFP»LtGVAR(P5l 11)1 
Tf(IERM.Nf.O)  go  to  <,00 
100  p*IMAx  a  «M«xl(An,(P^I(T)).PeIMA|0 
JSloP  ■  I  STOP  /  P 

c 

C  PFCROEW  psi 

c 

no  200  I  ■  1,  J«tqp 

TPOFXalSToP*!-! 

TrPpsPSl ( T 1 

os  1 (  1) *PS T  I  TN0F»' 

200  P«I  ( INUEX)«TFMP 

pfT|J(Jn 

sop  continue 

If  (  lEPrl  ,|  T,  o'  TFRR«S 

oaLL  EHRPoO(IE»o  ♦  Vi 

fnO 


» 


* 
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enonruilNE  tnTnp 


1  *n£CK  TKTNP 

el,HQ0ulINF  IKtMO 

rr»'M0N/  INPilTP/H.Mfj*  *o  «Fo5*FT0PmU»STqPk«  STOP* 
rr>*M0N  PSf  (50  01  .M (FO o  J  f «*0 «nr|.  Ta**NT.JTITLF 
5  rr^MON  /Fiino/  apf 

OpAL  NO.N.MIJ  IK 
pxTfRNAL  FUNcT?.  F'lNrT* 

/inG  ,  * 

CI JN 3  ■  .5»  t  (M  1 1  «PFl  li>  I  •*?  ♦  tN(NO)  «P$!  (noi  I  **2) 

io  tfPm»mU»r»mp 

T«rnP»No_i 
rn  100  I*?«!«Too 
100  ci)MsSuM3* <n (  I) •  o$f ( T ) ) *#2 

CiiKT  s  Su^T  •  nrl  x 4 x 
15  TFSf  *  <2.  *  K  •  B>  *•? 

TPST  ■  ABPMTE5T  -  >.)  *  <Tf«T  -  B.  )  > 

TF?T  *  TEfT  /  1PB,  *  1.0 
TfST  «  SQBT ( TEBT i 
i  ipPpR  ■  AmJNI  <  TfFT •  TEPu) 

20  7 3  gaUSSU  .  0..  UPPfB.  FUm*T2> 

IF  (TE°m  .GT,  TffTjj  a  7  *  0*US*>(*»  TEPT.  TtRM,  F,,NcTa) 

Z  *  1  '  Mu*«? 

FIIP3»sUm3*7/b 

PFTi'iBN 

?5  FNC 
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UM.IICN  rUNCT? 


FLECTION  FDNCT?/*) 

C('f'MCN/INP|lTP/n.MO,XS,Fr‘‘;*^TrPMu,^TOPx*  SfUP* 

COMMON  PSt^flOl'NtSOn^xO^PPLTiXtNT.lTITLF 

CPMON  /FllMr/K 

np«L  NO,N,k 

CALL  (K«8, i , x 

PI!NcT2“BE<;<:#*2«* 

RfTiibn 

FK  0 
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FUM.IIPN  ffUNC T 4 


function  fMNCT*<») 

/•CPmON' INPUTS'^*  NO  •xS‘P°P,,TcP*U»SToPk*  5 TOP* 
COMMON  Ps T  t^OOT  ,M(^0o»  t^O.riFI  T*X.Nt,ItItLf 
rr^MON  /FlINCT/  te 

pf»l  no»w,k 


jiFYmPTUT  Ic  4PPB0*  JMATIOO  NOB  J  (  K#R  ,  X(##2#* 


PT,  ACOSc-1., 

TFfiMl«<i.  /PJ 
TfBm2«K*0 

&BG«X  -  |,P  •  ?*BPB  ♦  .’M  *  PI 

TfBh3  ■  (*.  •  TrOM?  ••  5  “  i . >  /  (P.  •  X> 

FI.'NfT*  ■  Tf-Bol  •(C^SfABnl  •TeHM3  *  «5lN  I*Nfl>  1 

BftiiBN 

*>C 
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cnpPouiUE  tnitIal 


1 


in 


15 


?n 


?5 


10 


15 


<UH0CI)TINf  IMTTAL  (PtJMi) 

Cn»'M0N/INPUT5/B.K/O«XS*E  =  ,'»5TPHMI)*SToPK*  STOP* 

prPMON  PSt  (5(10)  .kj(50n)  ,MO.nFlTAX*fcT«ITlTLF 
OPAL  N0,N 

c 

c 

C  CALCULATE  nELTA*  FPOm  X5  ANn  NO 

c 

c 

OFLTaX  ■  *n  /  (Nn  -  P) 

c 

c  PFAn  NUMEMCAL  waU'IEf 

c 

T«TnP  ■  No  -  1. 

PFAn(5*45?l tTITl F 
*5n  FrPMATtA^> 

PFAO(5»40O)  {*(’».!  >  1.TSTOP, 

c 

C  CALCULATE  Xfl*M  uaL'IE  OF  N 

C 

NO  ■  N(N0-1« 

MAo>  *  NO  *  E*oi-nEl  Ta*Xh» 

C  CALCULATE  THE  INTERNAL  "F  n 

C 

FIIF1.  »5*(N(1  )  *M,Nn)  , 

TATnp.*  Nn  -  i 
no  i0n  I  a  2«  I«T0P 

300  FU*i  *  SUU1  ♦  W*T) 

Fl'Fj  «  SUmj  •  nn  TA* 

F(|M  ■  0  *  NO  ♦ 
nFTijpN 

c 

c  FPPOP  PR0CE55INP 

C 

700  CALL  E^RPoO  <  t ) 
pfTiibn 

*00  FoRMATlSFifl.B) 

FN0 
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( 


chpOOUIIkE  nuTPUTK 


•necK  output 

«ll6o0uUNP  OiitP|!tK  J tcOu^T > 

rrVMON'INPllT5/p.MD,X;»FO?*«TpPxU»SToPK«« InK* 

COMMON  PSf  (SflOl  . «g  f  >  .MOtriFL  WX«MT,  I  TITLE 
I*TfGEh  STOPnU 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

r 


pf*L  Not  N 

common  /Out/  Xt*rtO>»  Y<ooo> 
fP(KOUNT  Nf  1,  an  TO  PO 


MOITE  iOEwTJFVT^n  PAOAMfTEH*  To  TAPE  7 
FOLLOWED  RT  the  v*(.UfS  of  MIX) 


“PITE (2,96) P.Xn.oTnP*, EOF, Nn.STOPMii 
TFToP  *  No  -  l 

io|te<2«9f>  (N(»i »T«i*I«TOP) 
pfThPN 
5o  CONTINUE 


POINT  NUMqEP  Of  PATrF  A»'D  I.TST 
T me  Pa*? IS  TO  OUTPUT 

woite  tHE  nijmbfo  OF  PATHS  and 
TME  PA«S  to  TaPf  2 


POINT  98,kOUNT 
uo I IE  <  2«  99)  KOUNT 
no  10  I*l,EOUNT 
upJjE (2f lPfl)  X ( T )  *  y ( T ) 
IQ  PSlixT  lot  XlT^VIt) 

9R  F^mJV(|HY»  15) 

96  FOPm*T<*fTS.7*PtO//) 

9§  FORMAT  jBFjo^j 


101 


FOOMAT<lX,?f?2»T» 


.n  nm'&i'i 


EmO 


j 

i 


.1 

j 


! 
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subrcuiIkE  furpr* 


l  »* 


15 


2ft 


25 


SlJPROUTlNF  EPRPon  (N) 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


This  EhROB  PROcr«STNft  S'lRNOllTINE  PRESuPPOSts 
AK  ICA  TYop  pNV f oQnMfNT  WITm  CUC  6*60  EOUIH. 
"tMi  N0S/RF1  OpfRaTtNG  STATE* 


THE  EXTERNAL  RFrFRFNrE  aBhtjpb  tS  A  RouTlNt 
hr  I  T  T  Eft*  I  At  compass  T  o  ReNtR  a  Tt  FRROR  meSSaofSi 
IMTIa'E  TRACE  back*  ANpi  AHOPI  THE 
JOB  HjlH  PUT  A  fittHP 


1005 

'OOfe 


TO  l7Q«.8oo,onn,9()2,T001  , 1 002.1003,1 00* ♦ 1 00^,1 00^. 100T>n 
CALL  aBRijOrus  L  mOpe  tHak  *oo  *A  t  RS  hfhf  NfcEoEn  FOR  THIS 


Gft 


1  OOT 
700 
BO? 
900 
r02 
1001 

1  002 
1OO3 

100* 


CALL  aBRTjOR(6<’  L 


'00 

thf  m«* ihiih  k  value  ok  T«pts  is  less 
thfRE  ahe  not  enough  Hunts  on  tapes) 


136  t 


THIS  roaNCH) 
than  THE  k 


/  0  ) 


CALL  ABRT jOr ( 2 o i  ENO-OF-FILF.  UNIT  5) 

CALL  ABRT  iOR<16'  SOua'T  tOO  I  AHOE  ) 

CALL  ABRTjORI  ?«l  INpiEFTNITf  OPEPAN0  0 
t-ALL  AggTjOfll  )•»!  OIVISTON  BY  ZERO) 

ALL  ABRtJOB(3oi  CaNnOt  FINn  FOR  tHTS  houEj 

TOO  HanT  POINTS) 


CALL  AtJRT  JOB  <  161 

CONTINUE 

CALL  ABRt,|OR  (33t 

FK'O 


Kn  and  K1  ON  SAME  ?IOE  of  CurvE  ) 


SHB«fUU*>f  OLOTMPf 


1 

5 


1  0 


15 


?0 


25 


30 


erpooll* I*F  PLOT-no (NBTF.MOriF  *Xk  .  XOm.  STEP » 
C/'»<MCN/lNP|iTs/P.*>’0,xS.fo«;,,TrPMii,«TpP((,  bT^P* 
frKMON  PSt  (Snot  .  M  ■  0  o  J  *‘'0 *pfl T  ax  »tot 

nt^fN^lON  nct(*».  *X f 11 1 •  yv(ll)*  xpLnT»So**i»  YPiOTtSOO) 

FOL  I  VAUENrP  (PST  .  vPI.OT  ) 

PAT A  Qi •  o . »  i«  / 

PFL<*>  *  FTfP* 

FALL  PSCAl  E  (DEL  (*n 
PFL 1 J )  »  STPpX  X  In* 
pp  100  I»i«ll 
**(  I)  •  (T  -  1)  •  nEl.  <?> 
vvtn  *  (T  -  6)  •  PEl  (4< 

100  rrNTlNUE 

nr  ?66  I»v*  nPt« 

XPLOT(I)  ■  (I  -  i>  •  DEi Tax 
200  continue 

CALL  nAxIS(0.,n>.iqHDST  Str /«0«T (MFtEO) ,-l * ,V.5,<>0 . , yy, 1 . ,*Mr6.1 , 


1  «) 

CALL  PLOT <fl.»  5.*  -li 

CALL  SYMBOL  (.5.  4.5,  .it  PMKat  O.i  2) 


fall  nUMBfB 


(*P.  4.5,  ,li  »«i  0.»  4HF 4 * ? * 
(.5.  4.*  .I',  bHPWESAai  6,.  At 
<1  .2.4.,.1,XfP»0.,4HF7.5,7i 


41 


C4LL  nAX!s(0,«  A.t^MOEPTH  T*  *ETEBS»  15»  in**  0.*  **•  l.»iNFTt2 

1  .  n 


CALL  SyhBpl<2.,.b.»  .2.  llnupoE  NUpbEb,  o.» 
CALL  NUMBfB  (4.A.-S*.  .»«  PonEt  0,,  2HI4*  J) 
Fa Ll  LlNElXPLOT.  YPLpT.  NPTq,  1,  0.  8,  rtD 


CALL  P*-0T  (12, 
PFTuBN 


11) 


FAC 


suporoi  inf  olote*5 


i 


«5 


10 


1«! 


?n 


?? 


30 


3? 


•nECK  Pl.CT 

«('8P0uTlNP  PL0t»'o<nPts,“,KMAx> 

COMMON'OUT/XPLOT/snO'  •  vPLOT(soe) 

ri^FN<!lON  del  <4t  .X*(il)Tvr  Ij  |» 

PpAL  kMaX 

. -  OP  TO  t lo*  ?0»  30) M 

10  rnOTINUg 

0PL(U  «  ", 
opl<2>  ■  atopk  /  in* 
nrL  (3)  ■  o. 

PEL  (4)  a  KM4, 

C4LL.PsC*I.E(nEL**n 

no  100  1  a  1 •  it 

xx  <  t  >  ■  (T-l)  •  nEL (?) 

YV<t)  ■  *1-1)  •  nELU) 

108  fit 

CALL  0*XI«  (  g,  •  n.  .  i^Hk  (oAUTuK'S  /  MET'iR )  •  19*  i  u.  »  0.  * xx  *  1  •  *4ME6.4*6) 
CaLl  5YH0OL(?..in,..Pl»o4MfMTfcHNAL  w/Vf  uiSPfcRSInN  ReLatTOn'S . 0, • 
•  34) 

caul  symbol (4., o.s,.?i,ttiti F.o.tiO) 

CAUL  DAXl«?<p,*n. .?hOmE6a*-«.)0**9o* *TY*| *.4M»6.4.6) 

.  RrTURN 

20  CALL  LlNE(xPUOt.YPLO_, No.-, i.0,0, DEL) 

CAUL  P»-OT(o..O..-n 

PfTijPN 

36  continue 

CALL  STmbol(4*5..S*«!>1»p1MM»),IMUH  ERROR  In  k  I?»n**2)) 

Pal. 

CAUL  wMERrixOLn.YOLD',F) 

CALL  NUMBER (XOLo.TOUn*, oi *«MAX*o.*PHEio*3*ln> 


Rr Turn 

99  fcPMAT<lX,E2n.1f»t 

ENO 
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oinrcruiiKE  bscale 


j 


s 


in 


is 


SURP01»TINC  PSCAI  F  (SCaLE( 
TrfiMa  i  • 

ipiscaCe  .ge.  i*"  no  TO  ? 

1  J«<CAt|»TfRM 

IFIJ.G'.OI  nn  Tn  in 
TpfiM  ■  T£bm  •  l". 

nr  fO  1 

?  j»ScalE/Tfbm 

IM  J  ;LE.  lnlnr  TO  »0 

tFBu»tERK#10, 

nr  TO  2 

in  scale»(J*.S)/(i«.«tEom> 
pfTiiPN 

20  SCALE»<J*.S1»Teom/i0. 

BFThSN 

fnO 
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